[TOC]
# 兽医宏基因组检测技术分析流程（V1.0,2025.6.15）
    # 版本Version: V1.0,2025.6.15
    # 操作系统适用: Centos7、Centos8、Centos stream9/10及Ubuntu22.04+等linux系统。
    # 作者：中国农业科学院北京畜牧兽医研究所_动物生物安全与公共卫生防控科技创新团队，沈青春

# 一、数据整理
## 1.样品数据整理与查看

```
mkdir -p  train2025/testdata    
cd train2025
# cp /home/rcgh/train2025/testdata/*  testdata/  # 仅在培训期间用于测试数据获取
# rename 'sub_' '' *.fq.gz                       # 修改文件名(CentOS)
# rename s/sub_// *.fq.gz                        # 修改文件名(Ubuntu)
```

# 二、质控
## 1.样品数据质控
### 1.1 质控前的质量评估

```
mkdir -p qc/html
conda activate qualify
fastqc -t 4 testdata/*.fq.gz -o qc/html
cd qc/html
multiqc ./*
cd ../..
```

### 1.2 质控

```
# 软件选择，根据自己的喜好习惯即可，均可达到目的。
# fastp软件，华大数据推荐使用。
mkdir -p qc/fastp_clean
ls testdata | cut -d '.' -f 1 | sed 's/..$//'|sort -u | while read id ;
do 
fastp -i testdata/${id}_1.fq.gz -o qc/fastp_clean/${id}_clean_1.fq.gz -I testdata/${id}_2.fq.gz \
-O qc/fastp_clean/${id}_clean_2.fq.gz --thread 3 -f 13 -F 13
done                                                                            #逐个运行，速度很慢，可以考虑使用并发运行。
```

```
# trimmomatic软件，Illimina数据推荐使用。
mkdir -p qc/trim_clean
trimmomatic PE testdata/Hebei_sheep_lungD_1.fq.gz testdata/Hebei_sheep_lungD_2.fq.gz \
                              qc/trim_clean/Hebei_sheep_lungD_1_paired.fq.gz qc/trim_clean/Hebei_sheep_lungD_1_unpaired.fq.gz \
                              qc/trim_clean/Hebei_sheep_lungD_2_paired.fq.gz qc/trim_clean/Hebei_sheep_lungD_2_unpaired.fq.gz \
                              ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:True LEADING:3 TRAILING:3 MINLEN:36
# 单样品处理，对于多样品可参照fastp的多样品处理方法完成。
```

### 1.3 质控后的质量评估

```
mkdir -p qc/html2
# conda activate qualify
fastqc -t 4 qc/fastp_clean/*.fq.gz -o qc/html2
cd qc/html2/
multiqc ./*
cd ../..
conda deactivate
```

# 二、去宿主

```
conda activate kneaddata
# mamba install kneaddata bowtie2 trimmomatic fastqc multiqc
kneaddata --version                           # kneaddata v0.12.2
trimmomatic -version                          # 0.39
bowtie2 --version                             # bowtie2-align-s version 2.5.4
```

```
# 创建子目录
db=/home/rcgh/DB
# mkdir -p ${db}/kneaddata
# cd ${db}/kneaddata
# # 下载
# wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/016/699/485/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_genomic.fna.gz  #鸡 (chicken)的参考基因组
# wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/016/772/045/GCF_016772045.2_ARS-UI_Ramb_v3.0/GCF_016772045.2_ARS-UI_Ramb_v3.0_genomic.fna.gz                        #羊（Ovis）参考基因组
# mv GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_genomic.fna.gz broiler.GRCg7b9485.2.fna.gz
# mv GCF_016772045.2_ARS-UI_Ramb_v3.0_genomic.fna.gz Ovis2045.2.fna.gz
# # wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/048/773/025/GCF_048773025.1_ASM4877302v1/GCF_048773025.1_ASM4877302v1_genomic.fna.gz                              #骆驼 (Camelus bactrianus)参考基因组
# # mv GCF_048773025.1_ASM4877302v1_genomic.fna.gz Camelus73025.1.fna.gz
# # gunzip Camelus73025.1.fna.gz 
# # bowtie2-build -f Camelus73025.1.fna Camel/Camelus73025.1 --threads 8
# # 解压
# gunzip broiler.GRCg7b9485.2.fna.gz  Ovis2045.2.fna.gz
# # bowtiew建索引，输入文件，输出文件前缀，4线程2分
# time bowtie2-build -f broiler.GRCg7b9485.2.fna broiler.GRCg7b9485.2 --threads 8  #13min
# time bowtie2-build -f Ovis2045.2.fna sheep/Ovis2045.2 --threads 8
```

```
# 建议多样品并行，占用磁盘空间约为原始数据5x空间
mkdir -p temp/sheep
time kneaddata -i1 qc/fastp_clean/Hebei_sheep_lungD_clean_1.fq.gz -i2 qc/fastp_clean/Hebei_sheep_lungD_clean_2.fq.gz \
      -o temp/sheep --output-prefix Hebei_sheep_lungD \
      --bypass-trim --bypass-trf --reorder \
      --bowtie2-options '--very-sensitive --dovetail' \
      -db ${db}/kneaddata/sheep/Ovis2045.2 \
      --remove-intermediate-output -v -t 10  
```  
      
```
#批量方法   
mkdir -p temp/chicken
for id in Ch_breast_good_98_RNA Ch_breast_sick_97_RNA ; do \
time kneaddata -i1 qc/fastp_clean/${id}_clean_1.fq.gz -i2 qc/fastp_clean/${id}_clean_2.fq.gz \
      -o temp/chicken --output-prefix ${id} \
      --bypass-trim --bypass-trf --reorder \
      --bowtie2-options '--very-sensitive --dovetail' \
      -db ${db}/kneaddata/chicken/broiler.GRCg7b9485.2 \
      --remove-intermediate-output -v -t 10   
    done 
```
  
```
    # mkdir -p temp/camel   # 骆驼
    # kneaddata -i1 qc/fastp_clean/Alq_camel_46_clean_1.fq.gz -i2 qc/fastp_clean/Alq_camel_46_clean_2.fq.gz \
    #   -o temp/camel --output-prefix Alq_camel_46 \
    #   --bypass-trim --bypass-trf --reorder \
    #   --bowtie2-options '--very-sensitive --dovetail' \
    #   -db /home/rcgh/train/trains/genomes/Camel/Camelus73025.1 \
    #   --remove-intermediate-output -v -t 10 
mkdir hostsq_rm
cp  temp/*/*paired_?.fastq ./hostsq_rm/  # 将去完宿主的序列收集到一起。
rename 'paired_' '' hostsq_rm/*.fastq    # 重命名
ll -h hostsq_rm/   
conda deactivate
```
# 三、注释

```
# conda env list
# conda create -n  kraken2 -y
conda activate kraken2
# mamba install -y kraken2=2.1.2 -c bioconda   #最新版本2.1.3运行结果有问题。

kraken2 --version # 2.1.2
mkdir -p temp/kraken2
db=/home/rcgh/DB/meta/kraken2/k2_veterinary/k2_veterinary202506
# db=/home/rcgh/DB/meta/kraken2/k2_veterinary/k2_chicken202506
```

```
ls hostsq_rm/Hebei* | cut -d '.' -f 1| cut -d '/' -f 2 | sed 's/..$//'|sort -u | while read id ; do
      kraken2 --db ${db} \
      --paired hostsq_rm/${id}_1.fastq hostsq_rm/${id}_2.fastq \
      --threads 8 --use-names --report-zero-counts \
      --report temp/kraken2/${id}.report \
      --output temp/kraken2/${id}.output; done
```

```
k2_vet1=/home/rcgh/DB/meta/kraken2/k2_veterinary/k2_chicken202506      
ls hostsq_rm/Ch* | cut -d '.' -f 1| cut -d '/' -f 2 | sed 's/..$//'|sort -u | while read id ; do
      kraken2 --db ${k2_vet1} \
      --paired hostsq_rm/${id}_1.fastq hostsq_rm/${id}_2.fastq \
      --threads 8 --use-names --report-zero-counts \
      --report temp/kraken2/${id}.report \
      --output temp/kraken2/${id}.output; done      
conda deactivate
```

# 四、结果呈现

```
# conda create -n krona
conda activate krona
# mamba install -y krona bracken krakentools r-optparse
mkdir report_all
cp temp/kraken2/*.report report_all/
ls hostsq_rm | cut -d '.' -f 1 | sed 's/..$//'|sort -u > metadata.txt
cp metadata.txt report_all/
cd report_all
# cat metadata.txt | while read id ; do  sed -i '1d'  ${id}.report; done
cat metadata.txt | while read id ; do kreport2krona.py  -r ${id}.report  -o ${id}report.txt ; done

# sed -i '' '1d' input.txt

cat metadata.txt | while read id ; do  sed -i '1d'  ${id}report.txt ; done
# sed '/某个字符串/d' file.txt  # 去掉含有该字符串的行。
cat metadata.txt | while read id ; do  sed -i '/Chordata/d'  ${id}report.txt ; done
ktImportText *report.txt  -o train2025_reports.html

conda deactivate
cd ..
```

# 五、宏基因组数据的拆分

```
## 需求：针对宏基因组检测数据中的某个感兴趣物种，提取其测序read序列，用于验证或深入分析。
# 从kraken2结果中提取目标物种的序列
```
## 5.1 安装taxonkit（推荐）
```
conda env list
# conda create -n seqtools
conda activate seqtools
# mamba install -c bioconda taxonkit seqtk
```

## 5.2. 确定目标物种的taxid

```
# 根据提示下载ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz，解压缩到 ~/.taxonkit/下。
cd ~/.taxonkit
wget -c ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz     #大小66.8M
tar -zxvf taxdump.tar.gz
# 测试
echo "Escherichia coli" | taxonkit name2taxid -r
echo "562" | taxonkit lineage -r
```

```
# 确定目标物种编号
echo "Mycoplasma capricolum" | taxonkit name2taxid -r       # 用于提取“山羊支原体”的核酸序列,2095。
echo "Mycoplasma" | taxonkit name2taxid -r                  # 用于提取“支原体”的核酸序列,2093。
```

## 5.2 从Kraken2输出文件中提取目标Read ID

```
# 设Kraken2输出的out文件，文件总共五列。
# 第一列：C/U，C代表分类，U代表未分类。
# 第二列：测序数据read的 ID号 ，来自FASTA/FASTQ文件。
# 第三列：分类ID号，NCBI的taxon id号;如分类id号为0,说明没有分类。
# 第四列：测序read的长度。
# 第五列：以空格分隔的列表，指示序列中每个k-mer的LCA映射 （冒号之前taxid，0为没有数据）。

mkdir temp/cf_Myco ; cd temp/cf_Myco
less ../kraken2/Hebei_sheep_lungD.output |grep "Mycoplasma capricolum" > M.cap_ids.txt
cat M.cap_ids.txt |grep "2095"  | wc

# echo "Orthoherpesviridae" | taxonkit name2taxid  # Orthoherpesviridae      3044472

# 提取所有子taxid对应的Read ID
# awk 'NR==FNR{a[$1]; next} ($3 in a)' M.cap_ids.txt ../kraken2/Hebei_sheep_lungD.output | cut -f2 > target_read_ids.txt
cat M.cap_ids.txt | cut -f2 > target_read_ids.txt           # 取对应的第2列，即read的ID。
# 5.3 从原始FASTQ中提取序列

# 使用 seqtk 根据Read ID列表提取序列：
# 单端测序数据（单端二代或三代数据）
seqtk subseq input.fastq target_read_ids.txt > output.fastq

# 双端测序数据（R1和R2）
# 在Read ID原始文件中带有 /1 或 /2 后缀,需要先添加（非常重要）
# 生成R1和R2的Read ID列表
sed 's/$/\/1/' target_read_ids.txt > target_read_ids_R1.txt
sed 's/$/\/2/' target_read_ids.txt > target_read_ids_R2.txt
# 提取序列
seqtk subseq ../../hostsq_rm/Hebei_sheep_lungD_1.fastq target_read_ids_R1.txt > M.cap_output_R1.fastq
seqtk subseq ../../hostsq_rm/Hebei_sheep_lungD_2.fastq target_read_ids_R2.txt > M.cap_output_R2.fastq

conda deactivate
```

# 六、组装分析 Assemble
## 6.1 去宿主（比较耗时间和计算资源） 
`    # 已完成，见“二、 去宿主”`
## 6.2 组装

```
# 启动工作环境
# mamba create -y -n megahit
conda activate megahit
# mamba install  -y megahit spades quast cd-hit emboss salmon prodigal
### megahit安装后测试
# megahit -v                                # MEGAHIT v1.2.9
# metaspades.py -v                          # SPAdes genome assembler v4.2.0 [metaSPAdes mode]
# metaquast.py -v                           # QUAST v5.3.0 (MetaQUAST mode)
# cd-hit -v | grep version                  # CD-HIT v4.8.1
# embossversion                             # EMBOSS v6.6.0.0
# salmon -v                                 # salmon v1.10.3
### MEGAHIT组装Assembly
# 删除目录中的旧文件夹，否则megahit无法运行
# 组装，TB级数据需几天或更长，MEGAHIT v1.2.9
mkdir temp/megahit
```

```
# 单个样品
id=Hebei_sheep_lungD
megahit -t 8 \
-1  hostsq_rm//${id}_1.fastq \
-2  hostsq_rm//${id}_2.fastq \
-o temp/megahit/${id}mega
```        

```
# 多样品循环
cat metadata.txt | while read id ; do time megahit -t 8 \
-1  hostsq_rm//${id}_1.fastq \
-2  hostsq_rm//${id}_2.fastq \
-o temp/megahit/${id}mega/; \
done
```      

```
# cp /home/rcgh/train/trains/seqkit ./
./seqkit stat temp/megahit/*/final.contigs.fa > megahit_results.txt  #统计结果
cat megahit_results.txt
```

```
# 预览重叠群最前6行，前60列字符
head -n 6 temp/megahit/Hebei_sheep_lungDmega/final.contigs.fa | cut -c1-60

# 备份重要结果
mkdir -p result/megahit/
cat metadata.txt | while read id ; do cp temp/megahit/${id}mega/final.contigs.fa result/megahit/${id}_final.contigs.fa ;done
    
# 删除临时文件
cat metadata.txt | while read id ; do  /bin/rm -rf temp/megahit/${id}mega/intermediate_contigs ;done

### (可选)metaSPAdes精细组装
# 精细但使用内存和时间更多，15~65m
mkdir -p temp/metaspades
cat metadata.txt | while read id ; do  time metaspades.py -t 32 -m 100 \
-1 hostsq_rm//${id}_1.fastq  \
-2 hostsq_rm//${id}_2.fastq  \
-o temp/metaspades/${id}meta/;\
done
```

```
ls -sh temp/metaspades/*/contigs.fasta                                  #查看拼接结果。
./seqkit stat temp/metaspades/*/contigs.fasta > metaspades_results.txt  #统计结果
cat metaspades_results.txt
# 备份重要结果
mkdir -p result/metaspades/
cat metadata.txt | while read id ; do cp temp/metaspades/${id}meta/contigs.fasta result/metaspades/${id}_contigs.fasta ;done
# 删除临时文件
rm -rf temp/metaspades
# 注：metaSPAdes支持二、三代混合组装

### QUAST评估
# QUAST评估，生成report文本tsv/txt、网页html、PDF等格式报告
cat metadata.txt | while read id ; do quast.py result/megahit/${id}_final.contigs.fa \
-o result/megahit/${id}_quast -t 2;\
done

cat metadata.txt | while read id ; do quast.py result/metaspades/${id}_contigs.fasta \
-o result/metaspades/${id}_quast -t 2;\
done  

```
 
## 6.3 基因预测、去冗余和定

```
### metaProdigal基因预测Gene prediction
# 输入文件：组装的序列 result/megahit/final.contigs.fa
# 输出文件：prodigal预测的基因序列 temp/prodigal/gene.fa
mkdir -p temp/prodigal
# prodigal的meta模式预测基因，>和2>&1记录分析过程提示或报错信息至gene.log
cat metadata.txt | while read id ; do time prodigal -i result/megahit/${id}_final.contigs.fa \
-d temp/prodigal/${id}_gene.fa \
-o temp/prodigal/${id}_gene.gff ;\
-p meta -f gff > temp/prodigal/gene.log 2 > 1 ;\
done
tail temp/prodigal/gene.log
```

```
# 统计基因数量
./seqkit stat temp/prodigal/*gene.fa > prodigal_fa.txt    #统计结果
cat prodigal_fa.txt                                       # 统计完整基因数量，数据量大时可只用完整基因部分
    
cat metadata.txt | while read id ; do grep -c 'partial=00' temp/prodigal/${id}_gene.fa ;done
```    

```
# 提取完整基因(完整片段获得的基因全为完整，如成环的细菌基因组)
cat metadata.txt | while read id ; do grep 'partial=00' temp/prodigal/${id}_gene.fa | cut -f1 -d ' '| sed 's/>//' > temp/prodigal/${id}_full_length.id ;done
cat metadata.txt | while read id ; do ./seqkit grep -f temp/prodigal/${id}_full_length.id temp/prodigal/${id}_gene.fa > temp/prodigal/${id}_full_length.fa ;done
./seqkit stat temp/prodigal/*_full_length.fa
conda deactivate
```    

```
### cd-hit基因聚类/去冗余
# 输入文件：prodigal预测的基因序列 temp/prodigal/gene.fa
# 输出文件：去冗余后的基因和蛋白序列：result/NR/nucleotide.fa, result/NR/protein.fa
conda activate megahit
mkdir -p result/NR
cat metadata.txt | while read id ; do cd-hit-est -i temp/prodigal/${id}_gene.fa \
-o result/NR/${id}_nucleotide.fa \
-aS 0.9 -c 0.95 -G 0 -g 0 -T 0 -M 0; done
```

```
grep -c '>' result/NR/*_nucleotide.fa
# 翻译核酸为对应蛋白序列, --trim去除结尾的
conda deactivate
cat metadata.txt | while read id ; do ./seqkit translate --trim result/NR/${id}_nucleotide.fa \
> result/NR/${id}_protein.fa ; done

```
        
# 7.耐药基因预测
## 7.1 CARD耐药基因

```
# CARD在线分析平台：https://card.mcmaster.ca/ 
# 本地软件使用教程: https://github.com/arpcard/rgi
# 参考文献：http://doi.org/10.1093/nar/gkz935
# 结果说明：protein.json，在线可视化；protein.txt，注释基因列表
mkdir -p result/card
# conda create -y -n rgi6
conda activate rgi6
# mamba install -c bioconda rgi

db=/home/rcgh/DB/card
# 数据库部署
# mkdir -p ${db} && cd ${db}
# 下载最新版数据库，40.5M (2025-5-29, 4.0.1)
# wget -c https://card.mcmaster.ca/latest/data
# 解压后35M
# tar -xvf data ./card.json
    
# https://www.jianshu.com/p/c6a72a523f75
# wget -c https://card.mcmaster.ca/download/0/broadstreet-v3.3.0.tar.bz2
# tar -vxjf broadstreet-v3.3.0.tar.bz2
    
# 加载数据库
rgi load --local -i /home/rcgh/DB/card/card.json
rgi database --version --local

# rgi load --card_json /home/rcgh/DB/card/card.json
# 宏基因组分析扩展数据库和加载
# rgi card_annotation -i card.json
# mv card_database_v4.0.1.fasta card.fasta
# rgi load -i card.json --card_annotation card.fasta

# 启动rgi环境和记录版本
conda activate rgi6
rgi main -v # 6.0.4
    
# 简化蛋白ID
cut -f 1 -d ' ' result/NR/Hebei_sheep_lungD_protein.fa > temp/Hebei_sheep_lungD_protein.fa
# 这个错误忽略即可，不是报错，没有任何影响  grep: 写错误: 断开的管道
grep '>' result/NR/Hebei_sheep_lungD_protein.fa | head -n 3
grep '>' temp/Hebei_sheep_lungD_protein.fa | head -n 3

# 蛋白层面注释ARG
# rgi load -i $db/card/card.json --card_annotation $db/card/card.fasta
time rgi main -i temp/Hebei_sheep_lungD_protein.fa -t protein \
-n 9 -a DIAMOND --include_loose --clean \
-o result/card/Hebei_sheep_lungD_protein
head -n3 result/card/Hebei_sheep_lungD_protein.txt
######----结果解读---#######
# ORF_ID: Piromyces_2569,表明检测到的 ORF（开放阅读框）编号。
# Start 和 Stop: 起始和终止位置，这里未提供具体的值，可能因为不是从全基因组数据中提取的序列
# Orientation: 基因方向信息（未显示，可能默认正向）。
# Cut_Off: Strict,指定匹配阈值，这里是严格（strict）模式，表示该抗性基因的注释达到了严格预测的要求。
# Best_Hit_Bitscore 和 Pass_Bitscore:Best_Hit_Bitscore: 21.9422,比对得分，用于衡量该序列与参考抗性基因的相似程度。分值越高，匹配越可靠。
# Pass_Bitscore: 500,这是检测通过的最低阈值，与实际比对得分比较（这里表明比对得分低于阈值）。
# Best_Hit_ARO: OXA-1228,表明最佳匹配的抗性基因是 OXA-1228。这是一个 OXA 型 β-内酰胺酶，与 β-内酰胺类抗生素（如青霉素和碳青霉烯类）的抗性相关。
# ARO: 3008638
# Antibiotic Resistance Ontology (ARO) ID，用于在 CARD 数据库中唯一标识该抗性基因。
# Drug Class: carbapenem; penicillin beta-lactam,该基因赋予对 碳青霉烯类 和 青霉素类 β-内酰胺抗生素 的抗性。
# Resistance Mechanism: antibiotic inactivation,该基因通过“抗生素失活”机制赋予抗性，通常通过水解抗生素分子。
# AMR Gene Family: OXA beta-lactamase; OXA-62-like beta-lactamase,所属的抗性基因家族，包括 OXA β-内酰胺酶和与 OXA-62 类似的 β-内酰胺酶。
# Predicted_DNA 和 Predicted_Protein:Predicted_DNA 提供核苷酸序列。
# Predicted_Protein 提供翻译的氨基酸序列。
# CARD_Protein_Sequence:CARD 数据库中的参考蛋白序列，可用于进一步分析基因功能。
# Percentage Length of Reference Sequence: 18.49,该序列覆盖参考抗性基因的比例（18.49%），说明匹配区域较短。
# Nudged: TRUE,表示该基因注释结果通过 "Nudge" 策略进行了调整，尽管严格模式未完全满足，但匹配接近严格模式要求。
# Note: loose hit with at least 95 percent identity pushed strict表明该基因匹配结果属于宽松命中（loose hit），但因为匹配身份（identity）超过 95%，被调整为严格命中（strict hit）。
# Antibiotic: oxacillin,该基因与 oxacillin（苯唑西林） 抗性相关。这是 β-内酰胺类抗生素的一种。
conda deactivate
```
    
# 7.2、abricate耐药性分析
```
# 软件的搭建
# conda create -n antibioc
conda activate antibioc
# mamba install -y abricate                                  
abricate --list                                      # 查看现有数据库，BRicate可以调以下9个数据库的最新版本
# • ARG-ANNOT：基于argannot耐药数据库
# • CARD：基于加拿大CARD耐药数据库
# • EcOH：用于大肠杆菌O/H抗原分析
# • Ecoli_VF：基于大肠杆菌毒力基因
# • MEGARES：基于megares耐药数据库
# • NCBI：基于NCBI的耐药数据库 
# • Plasmidfinder：用于质粒扫描
# • resfinder：基于resfinder耐药数据库
# • VFDB：基于vfdb毒力基因数据库
mkdir -p result/abricat; cd result/abricat
mkdir -p plasmidfinder card vfdb
cd -
cat metadata.txt |while read id; do \
abricate --db plasmidfinder --mincov 80 --minid 90 --csv result/metaspades/${id}_contigs.fasta > result/abricat/plasmidfinder/${id}.plasmidfinder.result.csv>&1;\
abricate --db card --mincov 80 --minid 80 --csv --csv result/metaspades/${id}_contigs.fasta > result/abricat/card/${id}.card.result.csv>&1; \
abricate --db vfdb --mincov 80 --minid 80 --csv --csv result/metaspades/${id}_contigs.fasta > result/abricat/vfdb/${id}..vfdb.result.csv>&1;\
done
# --minid DNA identity阈值设置，默认是80%
# --mincov DNA coverage 阈值设置，默认是80%
   
# #打包处理结果,用于下载。
# cd result/abricat
# tar -czvf vfdb.tar.gz vfdb/
# tar -czvf card.tar.gz card/
# tar -czvf resfinder.tar.gz resfinder/ 

conda deactivate
```

#参考资料：
## 1.Yong-Xin Liu, Yuan Qin, Tong Chen, Meiping Lu, Xubo Qian, Xiaoxuan Guo, Yang Bai, A practical guide to amplicon and metagenomic analysis of microbiome data, Protein & Cell, Volume 12, Issue 5, May 2021, Pages 315–330, https://doi.org/10.1007/s13238-020-00724-8
## 2.Garfias-Gallegos D, Zirión-Martínez C, Bustos-Díaz ED, Arellano-Fernández TV, Lovaco-Flores JA, Espinosa-Jaime A, Avelar-Rivas JA, Sélem-Mójica N. Metagenomics Bioinformatic Pipeline. Methods Mol Biol. 2022;2512:153-179. doi: 10.1007/978-1-0716-2429-6_10. PMID: 35818005.
## 3.Lu J, Rincon N, Wood DE, Breitwieser FP, Pockrandt C, Langmead B, Salzberg SL, Steinegger M. Metagenome analysis using the Kraken software suite. Nat Protoc. 2022 Dec;17(12):2815-2839. doi: 10.1038/s41596-022-00738-y. Epub 2022 Sep 28. Erratum in: Nat Protoc. 2024 Aug 29. doi: 10.1038/s41596-024-01064-1. PMID: 36171387; PMCID: PMC9725748.
